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Abstract 

^yt^ [ Formation of clusters of interstitial carbon atoms between the surface and second atomic layers of 

graphite is demonstrated by means of molecular dynamics simulations. It is shown that interstitial clusters 
result in the dome-like surface features that may be associated with some of the hillocks observed by STM 
on the irradiated graphite surface. 
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g 1. Introduction 

, Ion irradiation of graphite surface results in formation of defects that are seen by scanning tunneling microscopy 
(STM) |lU|. A single -ion impact creates a hillock whose nature still remains controversial since STM probes the 
electronic states of the surface rather than the actual surface topography. So, it is difficult to deduce the atomic 
structure in the vicinity of defects based on STM data alone. 

A number of suggestions have been made concerning the interpretation of hillocks seen by STM on graphite surfaces 
irradiated by low-energy (~ 100 eV) ions, including, e.g., the incident ions trapped between the surface and second 
layer of atoms in graphite |^|||, the vacancies on the graphite surface J^,[l0|,[ll) , self-atom implantation (jflj, etc. It 
seems, however, that in many cases the formation of hillocks involves carbon atoms only. 
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— i , Recent molecular dynamics simulations |11| indicate that carbon interlayer interstitials below the surface layer 
account for w 50 % of all defects produced by incident and secondary recoils with the energy E — 300 eV -j- 3 keV, 
while « 100 % of defects produced by secondary recoils with the energy E = 100 eV -j- 300 eV are single carbon 



interstitials. It was noticed in [fLll that the single interstitials may migrate to form clusters between the layers, which 



, may be the source of some of the hillocks. The authors of 1 1 1 arrived at a conclusion that, at most, half of the surface 
• defects may be formed by interstitial carbon clusters. 
ON , As to high-energy (~ 100 keV) ion implantation, the hillock shape of the defects was attributed by several authors 
to the formation of single carbon interstitials and interstitial carbon clusters between the graphite layers, each hillock 
being originated from a single ion impact ^,0. Two types of hillocks were found in ^J7|, one of which showing 
the undisturbed (i.e. definitely without vacancies) lateral atomic arrangement on the (0001) graphite surface, while 
the other one exhibited a distorted atomic structure (but also presumably without vacancies). So, it is likely that 
the hillocks formed under both low- and high-energy ion irradiation of graphite are mainly due to interlayer carbon 
interstitials and/or their clusters. If this explanation is true, a unique opportunity appears to study the processes 
of cluster formation and a related phenomenon of radiation-induced swelling of graphite by means of STM. Besides, 
there is a need for a theory of formation, structure, and arrangement of such clusters. The theory should give an 
answer to several key questions: 1) what is the reason for attraction of interlayer interstitials; 2) what are the physical 
conditions of cluster formation; 3) why the recombination of interstitials with vacancies does not inhibit the formation 
5J] ■ of clusters. 

The purpose of this paper is to study numerically the formation of clusters of interstitial carbon atoms in the 
interlayer region between the surface and second graphite layers. Recently we have shown by means of molecular 
dynamics simulations that interlayer carbon interstitials do form clusters in the bulk of graphite ]1^| . The physical 



reason for this effect |13 is the deformation interaction of interstitials mediated by the nearby graphite sheets. An 
interstitial stretches the lattice in its vicinity, and the stretched domain appears to be attractive for other interstitials 
located within the same interlayer region, see Fig. 1 . In other words, the formation of interstitial clusters minimizes 
the deformation energy relatively to the case of spatially separated interstitials due to minimization of the overall 
buckling of adjacent graphite sheets. It is important that mobility of interstitial clusters is much lower than that of 
single interstitials, and hence the recombination of interstitials with vacancies is strongly suppressed [|14[. 
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On general grounds, one can expect the same mechanism of clusters formation to be realized near the graphite 
surface. Indeed, the formation of clusters of noble-gas atoms inserted between the topmost and second atomic layers 
of graphite has been numerically demonstrated by Marton et al. PJl5| . We are not aware of similar studies on the 
clusters of carbon atoms. 

The paper is organized as follows. First we present the general theory of diffusion instability of the uniform 
distribution of defects in solids. Next we numerically simulate the formation of interstitial carbon clusters underneath 
the surface layer of graphite. Our calculations show that such clusters can indeed be formed. Their specific shapes 
are dictated by the initial distribution of interstitials in the interlayer region and by the strong covalent interaction 
between the interstitials. The interstitial clusters give rise to dome-like features on the graphite surface that are 
similar to those observed experimentally by STM. 



2. Diffusion instability of a system of defects in solids 

In this Section we sketch the main features of the theory of diffusion instability that can lead to formation of clusters 
of defects in disordered solids, see Jl3],[l4| for more details. We shall demonstrate that the uniform distribution of 
defects is unstable with respect to separation in defect-rich and defect-poor regions if the concentration of defects 
is sufficiently high and exceeds some critical value that depends on material parameters and temperature. The 
mechanism of instability is rather general: the deformation of the crystal lattice by defects results in anomalous flow 
of defects along the gradient of defect concentration n(r,t), i.e. opposite in direction to the diffusion flow of defects. 

In the simplest case of defects of the same type, the kinetic equation for n(r, t) is 

— = Q - R(n) - divj, (I) 

where Q is the source of defects, the term R(n) accounts for recombination processes, and j is the defect current 
density. In the case of isotropic continuum media, the expression for j reads 

J ~ dr + 3T dr ' [ ' 

where the first term corresponds to the diffusion of defects, while the second one describes the motion of defects in 
the field of elastic stresses an = diva. Here D is the diffusion coefficient and Vl is the dilatation volume (f2 > for 
interstitials and Q, < for vacancies; |f2| ~ a 3 , where a is the lattice constant). 

In turn, average elastic stresses (an) caused by defects satisfy the following equation of the theory of elasticity: 

where 

K* = K 11 4- 

3K 

K and fx are elastic moduli. Substituting Eq.(|3|) in Eq.(|2|), we obtain the following expression for j: 

It follows from Eq. that the current caused by elastic stresses is directed opposite to the diffusion component of the 
current. This result has a simple physical meaning. Let us consider, for example, the system of defects with Q > 0. 
Such defects (e.g., interstitials) stretch the lattice, i.e. increase the average lattice constant, see Eq.([|). But since 
interstitials themselves tend to move to the stretched regions of the sample, there is an anomalous current along the 
gradient of defects concentration. This is what we call the deformation interaction of defects. 
The steady-state spatially- uniform solutions of Eqs. (0), (||) appear to be unstable at 

Hence, in the case that the average concentration of defects exceeds the critical concentration n c , the agglomerations 
of defects (clusters) are formed. For T = 300K and typical values of Vl w 5 • I0 -23 cm 3 and K* k, 10 12 erg/cm 3 we 
have n c ~ I0 19 cm , this value being much smaller than the atomic concentration. 
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We stress that the results presented above are obtained within the framework of a standard defect model in the 
context of the linear theory of elasticity, without use of any new assumptions and parameters. Note also that the 
diffusion instability is a common feature of a number of other physical systems, e.g., superconductors and excitonic 
insulators |ll|. 

Our results hold true also for defects with f2 < 0, e.g., vacancies. Indeed, a vacancy compresses the lattice, i.e. 
reduce the average lattice constant, while the compressed domain appears to be attractive for other vacancies. If 
both defects with > and < (e.g., interstitials and vacancies respectively) are present in the sample, their 
effective interaction with each other via the lattice deformation is repulsive Q. Such a repulsion leads to spatial 
separation of vacancies and interstitials. Interstitials form the interstitial clusters, while vacancies form the vacancion 
clusters. The interstitial clusters are much less mobile than a single interstitial. As a result, the recombination is 
suppressed, and the concentration of defects under the influence of an external source can increase monotonously for 
a very long period of time without saturation Il4j . We stress that our model was formulated for the case of a finite 
defect concentration, and hence does not imply that, e.g., a single vacancy and a single interstitial have a repulsive 
interaction. 

It should be noted that while in isotropic crystals a characteristic range of deformation interaction is very short, 
about one lattice constant (and hence such an interaction can be viewed as point-like in the continuum limit fl3| , |l4|| ) , 
in anisotropic crystals the range of deformation interaction can be much longer. In particular, interlayer interstitials in 
graphite (both impurities and carbon atoms) are known to deform atomic layers very strongly JI^JitJ . As a corollary, 
the deformation interaction between the defects in graphite should play an important role in temporal evolution of the 
system of defects under irradiation and influence the spatial arrangement of defects over the sample. The formation of 
spatially separated interstitial and vacancy clusters and the corresponding suppression of recombination can explain 
some peculiar features of graphite swelling under irradiation Jl3| , pT[ . Note that other mechanisms (e.g., the specific 
bonding characteristics of the defects) can also contribute to the recombination suppression in graphite. 

The effect of interstitial carbon atoms formation in the bulk of graphite has been demonstrated numerically in |l2] ] 
making use of empirical interatomic potentials. Below we study the formation of such clusters beneath the surface 
graphite layer taking into account the strong covalent bonding between carbon atoms constituting the cluster. 



3. Computational details 

Graphite was modeled by six layers of hexagonal carbon networks, each layer staggered with respect to its neighbors, 
see Fig. 2. The layers were composed of 220 atoms each. The atoms at the periphery of each layer were fixed in order to 
keep the stability of the crystallite. All other atoms were allowed to move after the interstitials have been inserted into 
the crystallite. Initially (t = 0) ten interstitial carbon atoms were positioned randomly between one of the peripheral 
layers (the "surface layer") and the penultimate layer. 

In the presence of interstitials, the undisturbed graphite layers correspond to a highly nonequilibrium atomic 
configuration. The potential energy of the system rapidly decreases in time, leading to an increase in the kinetic 
energy. An additional decrease of the potential energy occurs because of the formation of covalent bonds between the 
interstitial carbon atoms upon fusion of interstitials into clusters. In order to avoid the overheating of the crystallite, 
we eliminated the whole kinetic energy of atoms artificially each 70 MD steps (each period of time t a = 2 • 10 -14 
s), i.e. in fact we removed the heat from the crystallite. A rough estimate of the system temperature made on the 
basis of values of the maximum kinetic energy attained each period i , gives T ~ 1000 K at the very beginning of 
simulations and T ~ 10 K at the final stage (after the formation of a single interstitial cluster). 

We made use of empirical interatomic potentials proposed by Taji et al. ]l7| to account for interaction between the 
carbon atoms located in graphite layers as well as between the interstitial carbon atoms and the atoms of graphite 
layers. Parameters of those potentials were chosen by fitting the calculated values of interatomic distances and elastic 
moduli to their experimental values in defect-free graphite. There arc two parts in the Taji potential. The first part 
describes the interaction between the atoms within a graphite layer and takes the non-central force interaction into 
account in the form of the Keating strain energy potential. This part has an explicit minimum for angles between 
C-C bonds of 120° and hence treats sp 2 bonding. The second part describes the interaction between graphite layers 
by the central force potential. This part includes the attractive van der Waals force and the repulsive force that 
prevents graphite layers from merging together. The second part of the potential was applied by Taji et al. to the 
interactions not only between the carbon atoms located on the different layers but also between the interstitial carbon 
atom and other lattice atoms. For an interstitial located between graphite layers, its interaction with graphite layers, 
as described by the Taji potential, appears to be repulsive, i.e. the interstitial atom does not form covalent bonds with 
carbon atoms of the nearby graphite layers. Such a choice of the potential agrees with the experimentally observed 
expansion of graphite due to interstitials that gives evidence for the repulsion between interstitials and graphite 
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layers, i.e. for the absence of covalent bonds between the interstitial atom and the graphite layers. In other words, 
the interstitial self-energy Ei s is positive. Its value Ei S = 1.8 eV calculated using the Taji potential agrees well with 
experimental data and theoretical estimates [ ^8|Jl9| ] . 

Note that the Taji potential was originally proposed E^] for the case that the interatomic distances and bond angles 
are not so far from their equilibrium values in graphite. This condition is, however, fulfilled in the simulations of 
formation of both a single interstitial 1 17 and interstitial clusters since the distortion of the surface layer by interstitials 
is rather smooth, see below. 

Of course, the Taji potential is empirical and its use should be also justified by more sophisticated calculations. 
Since first principles calculations of interstitial characteristics in graphite are not known to us, we have carried out 
TBMD (tight binding molecular dynamics |2^,^l|) simulations of interstitials between the graphite layers. We have 
demonstrated that there is indeed a repulsive force between an interstitial and the nearby graphite layers, so that new 
covalent bonds are not formed and sp 3 bonding does not appear. We have also shown that there is not only qualitative, 
but also semi-quantitative agreement between TBMD calculations and simulations based on the Taji potential. In 
particular, the calculated values of the volume expansion due to an interstitial, the displacements of atoms in the 
nearest to an interstitial graphite layers, and the value of Ei S are close to those found using the Taji potential. 

In order to account for the covalent interaction between the interstitials we employed the TBMD method that had 
been proven to give a good description of small carbon clusters [ECN 03] . 



4. Results and discussion 

Fig. 3 shows the dynamics of 10-interstitial cluster formation in the interlayer region nearest to the graphite sur- 
face. Initially (t = 0) the distance between any two interstitials greatly exceeds the characteristic range of covalent 
interaction (« 2 A). At the first stage the motion of interstitials is governed entirely by the attractive deformation 
interaction fll2|] . At t = 100, a small 2-interstitial cluster is formed. The atoms of the cluster are tightly bound 
together by covalent bonding, the covalent component of the binding energy being Ef, — 2.85 eV/atom, in agreement 
with the experimental value for the dimer C2 fl23| . At t = 200, one more 2-interstitial cluster is formed. In fact, 
at t — 200 -7- 500 two processes take place, the formation of new diatomic clusters and adsorption of single inter- 
stitials by those clusters. The main reason for the motion of single interstitials towards the interstitial clusters is 
the deformation-induced attraction [Q, as in the case of 2-interstitial clusters formation. At t = 500, there are two 
3-interstitial clusters and one 4-interstitial cluster in the interlayer region. Those clusters have the linear chain struc- 
ture governed by the covalent interaction between the interstitials within the clusters pl| ]. The covalent component 
of the binding energy per interstitial is Et> = 4.73 eV/atom and Ej, = 4.75 eV/atom for 3- and 4-interstitial clusters 
respectively, in accordance with theoretical results of other authors and experimental values (see, e.g., ]24| , p5| ] ) . 

At the next stage of the evolution (t = 500 -j- 5000), three interstitial clusters slowly move towards each other. At 
t = 5000, two clusters merge into a 6-interstitial cluster. Note that characteristic times of 3- and 4-interstitial clusters 
formation are an order of magnitude shorter than the time it takes for those clusters to merge together. In other 
words, the mobility of an interstitial cluster decreases rapidly as the number of interstitials in the cluster increases. 

Finally, the 10-interstitial cluster is formed, see Fig.3h. The shape of this cluster is typical for low-dimensional 
carbon structures that are characterized by the bond angles 180° (carbyne) and 120° (graphite layers) and by coor- 
dination numbers Z = 2 and Z = 3. The covalent component of the binding energy of interstitials within the cluster 
is Eb = 5.50 eV/atom. Note that the most stable geometry of an isolated 10-atom carbon cluster is a ring structure 
|l|,|6|j27|, its binding energy calculated using the same TBMD method equals to E^ = 5.87 eV/atom. Nevertheless, 
we have verified that the shape and the binding energy of the 10-interstitial cluster shown in Fig.3h remained nearly 
unchanged upon its removal from the crystallite and subsequent relaxation. Hence this cluster, taken alone (outside 
the graphite) is metastable and corresponds to a local energy minimum. 

The deformation of the surface layer by the 10-interstitial cluster is shown in Fig. 4. One can see the dome-like 
feature on the graphite surface. Its height is about 1.5 A and lateral dimensions are about 10 12 A. We suggest that 
some of the hillocks seen by STM on irradiated graphite surfaces [|]-|| may be due to relatively large (composed of 
~ 10 carbon atoms) interstitial clusters formed in the interlayer region between the surface and penultimate graphite 
layers at the sacrifice of deformation interaction of interstitials. The case considered by us (10 interstitials are created 
simultaneously in the surface area ~ 10 A 2 ) obviously corresponds to rather high energy of incident ions. Note, 
however, that interstitials' dynamics shown in Fig. 3 may be viewed as the final stage of cluster formation from well 
isolated interstitials initially created far apart from each other by low-energy ions and migrated to the region confined 
by the crystallite due to deformation attraction and/or thermally activated diffusion. 

We have made molecular dynamics simulations of interstitial clusters formation for several other initial distributions 
of interstitials over the interlayer region nearest to the graphite surface. The results are shown in Figs. 5-7. They 
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are similar to those presented above. One can see that in either case studied all interstitials, initially located far 
apart from each other, merge into a single 10-interstitial cluster. The shape of the cluster is not unique, it obviously 
depends on the initial distribution of interstitials and hence on the specific sequence of small clusters formation and 
subsequent fusion into a single large cluster. However, the angles between the C-C bonds in those clusters are equal 
or close to either 180° or 120°, while the coordination numbers of interstitials within the clusters are Z — 2 or Z = 3 
(with obvious exception of boundary interstitials for which Z = 1), just as in the cluster shown in Fig.3h. Thus the 
shape of the interstitial cluster is strongly influenced by the covalent interactions between the interstitials within the 
cluster. 

The covalent components of the binding energy of interstitials with each other equal to Ef, — 5.34 eV/atom, 5.32 
eV/atom, and 5.68 eV/atom for the clusters shown in Figs. 5b, 6b, and 7b respectively. Those values of Ef, are lower 
than the binding energy Ef, = 5.87 eV/atom of a stable 10-atom ring. Nevertheless we have found that the clusters 
preserve their shapes upon their removal from the crystallite, the values of Ef, being changed insignificantly Thus all 
final configurations of interstitials are metastable. 

Distortions of the surface layer by 10-interstitial clusters are shown in Figs. 5c, 6c, and 7c. It is clearly seen that 
the shape of the dome-like feature on the graphite surface reflects the particular shape of interstitial cluster. The 
height of the hillock varies slightly from one interstitial configuration to another, its typical value being 1.5 A. The 
lateral dimensions of the hillock are more sensitive to the shape of the cluster, usually being in the range 10 -j- 15 A. 

Finally, we have also simulated the deformation of the graphite surface by a single carbon interstitial and by 
interstitial clusters composed of various numbers of interstitial carbon atoms. We have found that the calculated 
height of the surface hillock increases monotonously with the number of interstitials in the cluster, starting from 
about 0.8 A for a single interstitial. The same holds true for characteristic lateral dimensions of the hillock. We 
believe that such a variation in size of hillocks may be at least a part of the reason for experimentally observed variety 
of features on irradiated graphite surfaces. For example, Coratger et al. have found that irradiation of HOPG 
graphite with 20 keV 12 C + ions resulted in protrusions presented an elongated form. These finding can be explained 
by the formation of relatively small clusters Cat beneath the graphite surface since for N < 10 a linear chain structure 
of carbon clusters ia believed either to be the most stable geometry or to lie very close in energy to a ring structure 



5. Summary and conclusions 

Making use of molecular dynamics simulations, we have demonstrated that interstitial carbon atoms formed in the 
interlayer region between the surface and second graphite layers under ion irradiation attract each other. The physical 
origin for such an attraction is the strong deformation of graphite layers by interstitials that makes it energetically 
favorable for interstitials to come closer to each other (and thus to minimize the overall buckling of graphite layers) 
than to stay far apart. 

Attractive deformation interaction of interstitials results in formation of interstitial clusters from initially isolated 
single interstitials. The specific shape of the cluster depends on the initial distribution of interstitials over the 
interlayer region and is dictated by the strong covalent interaction between the interstitials. Coordination numbers 
for the majority of interstitials in the cluster are Z = 2 and Z = 3, while all angles between the C-C bonds in the cluster 
are close to 180° or 120°, these values being typical for low-dimensional carbon structures. Various arrangements of 
interstitials in the clusters correspond to different metastable states of the cluster at a fixed number of carbon atoms 
in it. 

Deformation of the surface layer by interstitial clusters results in the dome-like features on the graphite surface. 
The typical height and lateral dimensions of the hillocks produced by 10-interstitial clusters are 1.5 A and 10 -j- 15 
A respectively. Such features have been observed repeatedly by STM on the graphite surfaces irradiated with noble 
gas ions. The results presented in this paper provide an explanation to experimental observations. In order to make 
a closer comparison between the theory and experiment, it would be interesting to calculate STM images from the 
hillocks formed above interstitial clusters of different shapes in graphite. 
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FIGURE CAPTIONS 



Fig.l. Distortion of graphite lattice by two interstitials located between the nearest graphite layers at the in-plane 
distance r = 4.9 A from each other. 

Fig. 2. Side view (a) and top view (b) of the crystallite composed of 1320 carbon atoms (six layers 220 atoms each). 
Closed and open circles show the atoms of odd and even (from the top to the bottom) layers respectively. 

Fig. 3. Dynamics of 10-intcrstitial cluster formation in the interlayer region nearest to the graphite surface. Top 
view. Atoms of graphite layers are not shown. Time t is measured in units of t Q = 2 ■ 10~ 14 s. t = (a), t = 100 (b), 
t = 200 (c), t = 300 (d), t = 500 (c), t = 2000 (f), t = 5000 (g), t = 20000 (h). 

Fig. 4. Deformation of the surface layer by the 10-interstitial cluster. 

Fig. 5. Initial (a) and final (b) arrangement of 10 carbon interstitials in the interlayer region between the surface 
and penultimate layers of graphite. Top view. Atoms of graphite layers are not shown. Deformation of the surface 
layer by the 10-interstitial cluster (c). 

Fig. 6. Same as in Fig. 4 for other initial distribution of 10 interstitials over the interlayer region. 

Fig. 7. Same as in Fig. 4 for other initial distribution of 10 interstitials over the interlayer region. 
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